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Abstract 

A  steady-state,  three-dimensional  model  of  a  complete  polymer  electrolyte  membrane  fuel  cell  (PEMFC),  including  both  the  anode  and 
cathode,  is  formulated  and  solved  using  a  finite  volume  computational  fluid  dynamics  (CFD)  code,  Fuel3D,  developed  at  Foughborough 
University.  The  model  is  first  validated  against  data  obtained  from  the  literature  on  a  global  basis.  It  is  further  validated  on  a  local  basis  using 
experimental  data  obtained  from  a  segmented  cell.  Excellent  agreement  is  obtained.  The  validated  model  is  then  used  to  study  the  effect  of 
electro-osmotic  drag  and  diffusion  of  water  across  the  membrane.  Overall  transport  of  water  across  the  membrane  is  seen  to  take  place  from 
the  anode  to  the  cathode  side.  Finally,  the  model  has  been  used  to  carry  out  some  parametric  studies,  such  as  variation  of  electrode  thickness, 
shoulder  width,  degree  of  permeability  and  oxidant  concentration,  to  provide  a  clearer  understanding  on  how  changes  in  parameters  affect  the 
cell  performance. 

©  2004  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

It  is  known  that  water  content  plays  a  very  important  role  in 
determining  the  performance  of  fuel  cells.  Too  little  water  will 
lead  to  drying  of  the  membrane,  which  reduces  its  conduc¬ 
tivity.  Too  much  water  will  cause  flooding  in  the  electrodes 
which  hinders  the  amount  of  reactant  gases  from  reaching  the 
catalyst  sites.  In  order  to  understand  the  transport  phenom¬ 
ena  that  occur,  careful  modelling  of  the  fuel  cell  has  to  be 
done.  Recent  one-dimensional  models  of  the  PEMFC  have 
been  developed  by  Bernardi  et  al.  [1,2],  Springer  et  al.  [3], 
Eikerling  et  al.  [4],  Wohr  et  al.  [5],  Baschuk  et  al.  [6],  Rowe 
et  al.  [7]  and  Maggio  et  al.  [8].  These  models  are  useful  in 
providing  insight  and  reasonable  predictions  of  the  cell  per¬ 
formance  in  the  low  and  intermediate  current  density  ranges, 
but  fail  to  reproduce  the  abrupt  drop  observed  experimentally 
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at  high  current  density.  Two-dimensional  effects  are  found  to 
be  important  and  to  have  a  significant  impact  on  some  aspects 
of  fuel  cell  operation  and  water  management  in  particular. 
More  recent  two-dimensional  models  are  given  by  Nguyen 
et  al.  [9,10],  Hubertus  et  al.  [11],  Scott  et  al.  [12],  Squadrito 
et  al.  [13],  Gurau  et  al.  [14],  Yi  et  al.  [15],  Singh  et  al.  [16], 
Dannenberg  et  al.  [17],Hsing  etal.  [18],Umetal.  [19],Paola 
[20]  and  Siegel  et  al.  [21]. 

In  order  to  have  a  better  understanding  of  how  the  actual 
fuel  cell  performs,  it  is  necessary  to  have  a  model  which  is 
three-dimensional.  This  is  especially  so  in  the  electrode  since 
its  role  is  to  allow  a  spatial  distribution  of  the  current  density 
on  the  membrane  in  both  the  direction  of  bulk  flow  and  the 
direction  orthogonal  to  the  flow  but  parallel  to  the  membrane. 
It  is  also  important  to  include  the  anode  so  that  the  movement 
of  water  across  the  membrane  can  be  accounted  for. 

Recent  three-dimensional  works  have  been  published  by 
Shimpalee  et  al.  [22-25],  Berning  et  al.  [26-28],  Zhou  et 
al.  [29],  Jen  et  al.  [30],  Um  et  al.  [31],  Nguyen  et  al.  [32] 
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Nomenclature 

a 

water  activity 

A 

specific  surface  area  of  the  control  volume 
(m-1) 

Ci 

concentration  of  specie  i  (molm-3) 

£>h2o 

diffusion  coefficient  of  water  (m2  s-1) 

£>a,b 

binary  diffusion  coefficient  (m2  s-1) 

F 

Faraday  constant  (96,487  C  mol-1) 

I 

current  density  (Am-2) 

I0 

exchange  current  density  (Am-2) 

mt 

mass  fraction  of  species  i 

M 

mixture  molar  mass  (kg  mol-1) 

nd 

osmotic-drag  coefficient 

P 

pressure  (N  m-2) 

P 

water  saturation  pressure  (N  m-2) 

R 

universal  gas  constant  (8.314  Jmol-1  K-1) 

S 

source  terms  (kg  m-3  s-1) 

fin 

membrane  thickness  (m) 

T 

temperature  (K) 

u 

velocity  vector  (ms-1) 

V 

voltage  (V) 

V 

volume  (m3) 

mole  fraction  of  species  i 

X 

x  position  coordinate  (m) 

y 

y  position  coordinate  (m) 

z 

z  position  coordinate  (m) 

Greek  symbols 

a 

net  water  transport  coefficient  per  proton 

G 

porosity 

0 

overpotential  (V) 

K 

permeability  (m2) 

[l 

viscosity  (kgm-1  s-1) 

(P 

any  variable 

P 

density  (kgm-3) 

cr 

membrane  conductivity  (S3-1) 

Subscripts  and  superscripts 

a 

anode 

A 

species  A 

B 

mixture  B 

c 

cathode 

cr 

critical 

eff 

effective 

g 

gas 

h2o 

water 

i 

species 

K 

anode  or  cathode 

m 

membrane 

np 

non-porous 

n2 

nitrogen 

oc 

open  circuit 

o2 

oxygen 

p  porous 

sat  saturation 

Miscellaneous  symbols 
()  superficial  average 

( )g  intrinsic  average 


and  Senn  et  al.  [33].  These  models  have  only  been  validated 
with  experimental  data  from  global  polarization  curves  and 
as  such  are  not  validated  on  a  local  level,  i.e.  how  well  they 
predict  the  local  current  density  distribution.  In  fact,  most 
of  the  models  found  in  the  literature  are  shown  to  predict 
the  global  fuel  cell  performance  rather  well,  be  it  for  one- 
[1-8],  two-  [21]  or  three-dimensional  [22,23]  models.  This 
agreement,  irrespective  of  the  dimensionality  of  the  model, 
might  be  due  to  running  the  experiments  under  conditions 
that  correspond  with  the  model  assumptions/simplifications. 
However,  it  is  also  quite  likely  that  the  agreement  might  be 
due  to  parameter  adaption  to  experimental  data,  usually  in 
the  form  of  global  polarization  curves. 

In  view  of  the  number  of  unknown  parameters  that  arise 
in  fuel  cell  modeling,  the  approach  of  global  validation  of 
a  model  is  not  sufficient.  To  extend  the  work  on  three- 
dimensional  modeling  of  polymer  electrolyte  membrane  fuel 
cells,  a  model  for  the  whole  fuel  cell  is  herein  presented  and 
validated  against  local  current  density  obtained  experimen¬ 
tally  with  a  segmented  cell,  equipped  with  a  flow  field  com¬ 
prising  of  parallel  flow  channels.  Such  a  comparison  of  model 
predictions  with  local  experimental  current  density  lends  the 
model  a  higher  credibility. 

The  model  is  similar  to  that  of  Shimpalee  et  al.  [22-25], 
in  that  conservation  of  mass,  momentum  and  species  in  the 
gas  phase  are  considered  on  the  cathode  and  anode  side.  A 
membrane  model  derived  by  Springer  et  al.  [34],  together 
with  appropriate  boundary  and  constitutive  relations  close  the 
model.  Contrary  to  Shimpalee  et  al.  [22-25],  however,  only 
the  minimum  number  of  species  equations  are  solved,  thus 
reducing  the  computational  cost.  Furthermore,  the  porous  na¬ 
ture  of  the  electrodes  is  treated  with  proper  volume  averaged 
quantities  and  care  is  given  to  proper  coupling  of  these  with 
the  free  flow  in  the  flow  channels. 

After  a  short  introduction  to  the  governing  equations, 
boundary  conditions  and  closure  relations,  a  first  global  val¬ 
idation  of  the  model  implemented  in  the  in-house  code, 
Fuel3D  [35]  is  carried  out  by  comparing  a  simple  fuel  cell 
geometry  with  computational  and  experimental  results  from 
Shimpalee  et  al.  [22],  who  used  a  commercial  software,  Flu¬ 
ent.  Finally,  a  more  complex  flow  geometry,  taking  into  ac¬ 
count  local  validation  with  experimental  data  obtained  from 
a  segmented  cell  [36],  is  considered.  The  model  will  then  be 
used  to  give  an  insight  into  the  physical  behaviour,  the  dis¬ 
tribution  of  water  vapour  and  the  overall  cell  performance 
of  the  fuel  cell.  Further  parametric  studies  will  follow  based 
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on  a  simplified  geometry,  which  is  typical  of  a  geometrically 
repeating  domain  unit  away  from  the  end  effects  in  a  generic 
fuel  cell. 


2.  Mathematical  formulation 


The  assumptions  used  in  the  3D  model  are: 

•  steady  state  and  laminar  flow  of  gases; 

•  an  isothermal  cell  due  to  the  high  thermal  conductivity  of 
the  solid  materials; 

•  an  ideal  gas  mixture; 

•  the  catalyst  layer  is  homogeneous; 

•  presence  of  water  is  only  in  the  vapour  phase.  As  such, 
all  reference  to  water  in  the  following  text  refers  to  water 
vapour,  except  those  in  the  membrane  which  is  in  liquid 
form. 


The  governing  equations  consist  of  the  continuity,  mo¬ 
mentum  and  the  various  species  transport  equations.  For 
porous  regions,  superficial  and  intrinsic  properties  have  to 
be  introduced.  Superficial  averages  are  defined  as: 


(pdV 

v 


and  intrinsic  as, 


(1) 

(2) 


where  V  is  the  total  volume  of  the  representative  elementary 
volume  (REV)  and  is  the  volume  of  the  free  pores  in  the 
REV.  The  two  averages  are  related  through, 

(<f>)  =  e  (4>)(g)  (3) 


where  e  =  V^/  V  is  the  porosity.  To  save  on  notation  in  the 
forthcoming  model  description,  the  brackets  and  g  superscript 
on  intrinsic  properties  are  omitted. 


In  this  paper,  velocity  vectors  are  treated  in  the  superfi¬ 
cial  averaged  form  while  the  pressure,  density  and  species 
mass  fractions  are  treated  in  the  intrinsic  averaged  form.  The 
governing  equations  are  presented  in  Table  1 .  The  equations 
represent  the  case  of  flow  in  porous  media  when  0  <  e  <  1 . 
When  g  =  1,  they  represent  the  case  of  flow  in  non-porous 
media.  Four  species  are  considered,  namely,  hydrogen  and 
water  on  the  anode  side,  and  oxygen,  nitrogen  and  water  on 
the  cathode  side.  The  source  terms  Su2  and  So2  refer  to  the 
depletion  of  hydrogen  and  oxygen  at  the  anode  and  cathode 
control  volumes  next  to  the  membrane  respectively.  Sn2oa 
and  Sn2oc  take  into  consideration  the  production  of  water 
at  the  cathode  catalyst  layer  and  the  transportation  of  water 
across  the  membrane. 

Since  there  are  two  species  on  the  anode  side  and  three 
on  the  cathode  side,  it  is  only  necessary  to  solve  one  species 
transport  equation  on  the  anode  side  and  two  on  the  cathode 
side.  Here,  the  water  transport  equation  is  solved  on  the  anode 
side  while  both  the  oxygen  and  water  transport  equations  are 
solved  on  the  cathode  side.  The  water  transport  equation  is 
chosen  on  the  anode  side  because  it  gives  a  direct  coupling  of 
the  presence  of  water  between  the  anode  and  cathode  sides. 
The  other  species  can  be  obtained  from  the  constraint  that 
mass  fractions  must  sum  to  unity.  Hence, 

mn2  =  1  —  mH2oa  on  the  anode  side  (4) 

wn2  =  1  —  mo2  —  ^h2oc  on  the  cathode  side  (5) 

The  following  empirical  equations  are  used  to  calculate 
various  terms  in  the  electrochemical  model.  These  equations 
are  based  on  the  assumption  of  a  hydrated  Nation  117  mem¬ 
brane  and  are  taken  directly  from  the  work  of  Springer  et  al. 
[34].  The  net  water  transport  coefficient,  a,  is  given  by, 


a(x ,  y) 

=  nd(x,  y ) 


FDh2q(x,  y)[cH2ocCv  y)  CH2oaCv  y)] 

l(x,  y)tm 


Table  1 


Governing  equations 


Governing  equations 

Mathematical  expressions 

Source  terms 

Conservation  of  mass 

V  •  p(u)  =  Sca  or  Scc 

Sca  =  Srt  +  VtoOa  at  an°de/membrane  interface 

Scc  =  Sq2  +  ShoOc  at  cathode/membrane  interface 

SCa  =S cc  =  0  otherwise 

Momentum  equations 

72V  •  p(u)(u )  =  -5/p  +  V/xV (u)  +  Sm 

Sm  =  —  —K —  in  porous  region 

Sm  =  0  otherwise 

Hydrogen  transport  equation 

V  •  p{u)mn 2  =  VpD^ffBVmH2  +  Sh2 

Sho  =  ~  Mn15py  '>A  at  ano(je/membrane  interface 

Sht  =  0  otherwise 

Anode  water  transport  equation 

V  •  p(w)mH2oa  —  VpZ)fBVmH2oa  +  ^H2oa 

ce(x,y)MH  QI(x,y)A  . 

lSHo0a  = - f - at  anode/membrane  interface 

SH2Qa  =  0  otherwise 

Oxygen  transport  equation 

V  •  p{u)mo2  =  VpD^ffBVmo2  +  So2 

S0,  =  —  M°2^-y)A  at  cathode/membrane  interface 

So,  =  0  otherwise 

Cathode  water  transport  equation 

V  •  p(u)mu2Oc  =  V P^a,b  V mH2Oc  +  ^h2oc 

[1+2 a(x,y)\Mn  0I(x,y)A  .  ~ 

Sh0Oc  =  - Jf  - at  cathode/membrane  interface 

ShoOc  =  0  otherwise 
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where  (x,  y )  here  represents  the  position  on  the  membrane, 
Dh2o(x,  y )  represents  the  diffusion  coefficient  of  water  and 
cu2od(x,  y)  and  cu2oc(x > y)  represent  the  molar  concentration 
of  water  at  the  anode  and  cathode,  respectively.  I(x,y)  is  the 
local  current  density,  tm  is  the  membrane  thickness  and  F  is 
Faraday’s  constant.  The  first  term  on  the  right  hand  side,  n d, 
is  the  electro-osmotic  drag  coefficient,  describing  the  amount 
of  water  dragged  across  the  membrane  by  each  proton  from 
the  anode  to  the  cathode  side.  The  electro-osmotic  drag  itself 
is  a  function  of  the  activity  of  water  on  the  anode  side,  aa,  of 
the  MEA  (membrane  electrode  assembly,  which  consists  of 
the  anode  electrode,  anode  catalyst  layer,  membrane,  cathode 
catalyst  layer  and  cathode  electrode).  The  electro-osmotic 
drag  coefficient  is  calculated  as, 

nd(x,  y)  =  0.0049  +  2.02 aa  —  4.53^  +  4.09 a\  if  aa  <  1 

(7) 


n&(x,  y)  =  1.59  +  0.1 59(aa  —  1)  if  aa  >  1  (8) 

It  is  observed  that  at  high  current  density,  the  amount  of 
water  dragged  across  the  membrane  from  anode  to  cathode  is 
greater  than  the  back  diffusion  from  cathode  to  anode  (second 
term  on  the  right  in  Eq.  (6)).  This  results  in  a  net  transport 
of  water  from  anode  to  cathode  side,  resulting  in  possible 
dehydration  on  the  anode  side.  Hence,  it  is  reasonable  to 
assume  that  the  water  content  in  the  membrane  is  more  likely 
to  be  lower  on  the  anode  side  than  on  the  cathode  side,  and 
consequently  the  activity  of  water  on  the  anode  side  can  be 
used  to  calculate  the  electro-osmotic  coefficient  across  the 
whole  membrane.  The  water  diffusion  coefficient,  Dh7o  is 
calculated  from, 

£>h2o(+  v)  =  5.5  x  10-n«d  exp[2416«1/303)-(i/r))]  (9) 

The  water  concentration  on  the  anode  and  cathode  sides, 
cu2oa  and  ca2oc  in  Eq.  (6),  can  be  obtained  from, 

Ch2o,k(+  y)  =  2y^L(0M3  +  17.8aK  -  39.84 

^m,dry 

+  36.0<2j^)  if  <  1  (10) 

Ch2o,k(-x,  y)  —  (14  +  1.4(<2k  -  1))  if  0K  >  1 

^m,dry 

(ID 


(K=  either  a  or  c  for  anode/cathode  side,  respectively). 
Pm, dry  is  the  dry  PEM  material  density  while  Mm^ry  is  the 
equivalent  weight  of  a  dry  PEM.  The  activity  of  water  is 
defined  as, 


clk(x,  y) 


*h2o,k(*,  y)p(x,  y) 

psat 

rH20,K 


(12) 


where  p  is  the  cell  pressure  and  vh2o,k  is  the  mole  fraction 
of  water  on  either  the  anode  or  cathode  side.  The  saturated 
vapour  pressure  of  water,  which  is  dependent  on  the  temper¬ 
ature,  is  estimated  from, 

PfO  K  =  [0.00644367  +  0.000213948(7  -  273.0) 

+  3.43293  x  10"5(7  -  273.0)2  -  2.70381 
x  10“7(T  -  273. 0)3  +  8.77696  x  10“9 
x(T  -  273.0)4  -  3.14035  x  10_13(7  -  273.0)5 
+  3.82148  x  10_14(7  -  273.0)6]  x  1.013  x  105 

(13) 


Eqs.  (6)— (13)  allow  the  source  terms  in  the  species  trans¬ 
port  equations  to  be  calculated  (assuming  the  local  current 
density  is  known).  For  binary  gas  mixtures  at  low  pressure, 
Da,b  is  obtained  from  [37], 

_ pDa,b _ 

(7crA7crB)1/3(7’crA7,crB)5/12((l/MA)  +  (1/MB))1/2 

T 

y/  (  7c  r  A  TcrB  ) 

where  pc r  is  the  critical  pressure  and  Tcr  is  the  critical  tem¬ 
perature.  T  is  the  cell  temperature  and  M/  is  the  molar  mass  of 
the  species  considered.  The  values  of  a  and  b  are  taken  to  be 
3.64  x  10  4  and  2.334,  respectively.  The  effective  diffusion 
coefficient  applicable  to  the  porous  regions  will  be  given  by 
applying  the  Bruggeman  relationship  [38], 

<B=e%B  (15) 

The  mixture  viscosity  can  be  obtained  from, 


M  =  V.  , M4+  (16) 

where  the  individual  species  viscosity,  /x;,  can  be  obtained 
from  gas  tables  at  a  given  temperature.  In  this  work,  the  mix¬ 
ture  viscosity  is  assumed  to  be  constant  and  independent  of 
temperature  variation.  The  mixture  density  is  calculated  from 
the  ideal  gas  law, 


pM 

~RT 


(17) 


where  R  is  the  gas  constant  and  M  is  the  mixture  molar  mass 
which  is  related  to  the  individual  species  molar  mass  by  the 
following  equation, 


(18) 


The  species  molar  fraction  can  be  determined  from, 


M 


Xi  =  mi 


(19) 
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The  local  current  density  of  the  fuel  cell  is  given  by  the 
following  equation, 

I(x,  y)  =  ^(X-y)[V0C  -  Veen  -  !?(*,  V)]  (20) 

fin 

where  Voc  is  the  open  circuit  cell  potential,  Vceii  is  the  cell 
potential  and  om  is  the  membrane  conductivity,  which  is  a 
function  of  water  content,  and  is  determined  from, 


tfmO,  y) 


=  100 


0.00514 


M 


m,dry 


exp 


Pm,  dry 
[1268((l/303)-(l/D)] 


CH2oa(-y  y)  -  0.00326 


(21) 


The  cell  overpotential,  r\,  is  a  function  of  local  current 
density,  the  exchange  current  density  at  one  atmosphere  of 
oxygen,  and  the  partial  pressure  of  oxygen.  Only  the  cath¬ 
ode  overpotential  is  being  considered  here  since  the  anode 
overpotential  is  usually  an  order  of  magnitude  smaller  than 
that  of  the  cathode,  and  hence  can  be  neglected.  The  cell 
overpotential  can  be  calculated  from, 


RT 

ri(x ,  y)  =  -  In 

/v  'y)  0.5  F 


1.013  x  105/(jt,  y) 
IoPo2(x,  y ) 


(22) 


where  I0  is  the  exchange  current  density  at  one  atmosphere 
of  oxygen  and  Pq2  is  the  partial  pressure  of  oxygen  on  the 
cathode  side. 

In  the  current  CFD  code,  the  velocity-pressure  solution 
algorithm  used  is  the  steady  state  semi-implicit  method  for 
pressure-linked  equations  (SIMPLE)  algorithm  [39].  This  al¬ 
gorithm  is  essentially  a  guess-and-correct  procedure  for  the 


calculation  of  pressure  on  the  co-located  grid  arrangement. 
Since  the  code  uses  colocated  grids,  the  Rhie  and  Chow  [40] 
interpolation  method  is  used  to  avoid  the  problem  of  oscil¬ 
lations.  A  24  processor  Origin  2000  computer  is  used  to  run 
the  simulations.  The  solution  is  considered  to  be  converged 
when  the  relative  error  in  each  field  is  less  than  10~6. 


3.  Global  validation 

A  schematic  diagram  of  a  typical  fuel  cell  is  shown  in 
Fig.  1  and  the  geometrically  repeating  domain  unit  which 
can  be  identified  away  from  the  end  effects  regions  is  pre¬ 
sented  in  Fig.  2.  The  computational  domain  consists  of  an 
anode  gas  channel,  an  anode  electrode,  a  cathode  electrode 
and  a  cathode  gas  channel.  The  membrane  and  catalyst  layers 
are  so  thin  compared  to  the  electrodes  that  they  are  repre¬ 
sented  as  a  boundary  between  the  anode  and  cathode  porous 
electrodes.  The  anode/cathode  electrode  interface  is  a  bound¬ 
ary  on  which  conditions  representing  the  membrane/catalyst 
layers  will  be  applied.  The  boundary  conditions  are  given  in 
Appendix  A. 

No  distinction  was  made  between  intrinsic  and  superficial 
variables  in  [22] .  The  porosity  was  taken  to  be  equal  to  one 
and  only  the  diffusion  coefficient  of  each  species  was  reduced 
by  a  factor  of  0.5  to  account  for  porosity  and  tortuosity.  The 
same  will  be  applied  here,  but  only  to  allow  direct  comparison 
with  [22] .  In  general,  as  has  been  shown  in  [35] ,  it  is  important 
to  take  proper  account  of  the  nature  of  the  volume  averaged 
variables  being  used. 

For  this  simulation,  the  grids  were  fixed  at  10  cells  in  the  i- 
direction  (v),  10  cells  in  the  j-direction  (y)  and  200  cells  in  the 
k-direction  (z)  for  the  gas  channels,  while  for  the  electrodes, 


Coflow:  both  hydrated  fuel  (anode)  and  hydrated  oxidant  (cathode)  out 
Counterflow:  hydrated  fuel  (anode)  out  and  hydrated  oxidant  (cathode)  in  or  vice  versa 


Counterflow:  hydrated  fuel  (anode)  in  and  hydrated  oxidant  (cathode) 
out  or  vice  versa 


Top  view  of  a  single  fuel  cell 


Catalyst  layers  Gas  channel 


Front  view  of  a  single  fuel  cell 


Fig.  1.  Schematic  diagram  of  a  single  fuel  cell  configuration. 
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Fig.  2.  Schematic  diagram  for  a  repeating  unit  of  a  complete  anode/cathode  assembly  for  a  PEMFC. 


they  were  fixed  at  38  cells  in  the  /-direction  (x),  10  cells  in 
the  y-direction  (y)  and  200  cells  in  the  ^-direction  (z),  after  a 
grid  dependence  test  was  carried  out. 

In  order  to  compare  the  numerical  results  with  [22],  four 
different  cases  were  considered:  very  low  humidity,  low  hu¬ 
midity,  high  humidity  and  very  high  humidity  (on  the  cathode 
side)  while  all  other  parameters  were  kept  constant. 

Fig.  3  shows  the  comparison  of  results  between  those  re¬ 
ported  in  [22]  and  the  present  model.  The  operating  condi¬ 
tions  are  given  in  Table  2.  All  the  four  cases  considered  are  in 
coflow.  Current  density  is  presented  as  a  ratio  because  in  the 
work  of  [22],  the  experimental  data  were  obtained  based  on 
a  pressure  of  2  atm  while  numerical  predictions  were  carried 
out  based  on  an  operating  pressure  of  1  atm.  It  can  be  seen 


<u 


*  Degree  of  humidity,  1  -verylow,  2-low,  3-high,  4-veryhigh 

Fig.  3.  Average  current  density  ratio  dependence  on  the  inlet  humidity. 


that  the  results  obtained  from  the  present  model  are  in  ex¬ 
cellent  agreement  with  the  numerical  ones  from  [22]  for  all 
the  cases  considered.  Good  agreement  is  also  obtained  with 
the  experimental  results  obtained  in  [22],  except  for  the  case 
of  very  high  humidity.  As  the  degree  of  humidity  increases, 
the  average  current  density  increases.  This  is  due  to  the  fact 
that  the  amount  of  current  obtained  is  strongly  related  to  the 
degree  of  hydration  of  the  membrane  (Eqs.  (20)  and  (21)). 
As  a  result,  as  the  humidity  increases,  the  membrane  is  more 
hydrated,  resulting  in  better  performances.  However,  this  is 
only  true  if  the  electrode  is  not  flooded  with  excess  liquid  wa¬ 
ter  which  acts  as  a  barrier  to  the  transport  of  reactant  gases  to 
the  catalyst  layer  where  electrochemical  reaction  takes  place. 
This  can  be  seen  to  occur  in  the  experimental  result  at  very 
high  humidity.  Hence,  the  drop  in  current  density  from  the 
experimental  data  is  not  observed  in  the  simulated  result  since 
in  the  present  model,  water  is  only  considered  in  the  vapour 
phase. 

Another  point  to  note  is  that  in  [22],  five  transport  equa¬ 
tions  were  solved,  namely,  hydrogen  and  water  on  the  anode 
side,  and  oxygen,  nitrogen  and  water  on  the  cathode  side. 
However,  in  the  present  model,  only  water  on  the  anode  side, 
and  oxygen  and  water  on  the  cathode  side  were  solved.  The 
remaining  species,  such  as  hydrogen  on  the  anode  side  and 
nitrogen  on  the  cathode  side  can  be  obtained  by  the  mass  bal¬ 
ance  of  total  species  present  on  each  side.  This  will  be  more 
efficient  as  it  reduces  the  number  of  equations  to  be  solved 
numerically,  as  was  also  done  by  [19,26,29]. 
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Table  2 


Values  of  geometry  and  empirical  parameters  used  in  the  base  case  (reported 
data  used  in  [22]) 

Channel  width 

8.0  x  10“4  m 

Channel  height 

1.0  x  10-3  m 

Channel  length 

0.1  m 

Electrode  width 

3.2  X  10-3  m 

Electrode  height 

2.5  x  10~4  m 

Electrode  length 

0.1  m 

Membrane  thickness  (fm) 

5.0  x  10-5  m 

Permeability  of  electrode  (k) 

2.0  X  10-10  m2 

Porosity  of  electrode  (e) 

1.0 

Flow  conditions  at: 

Very  low  humidity 

Anode 

Inlet  velocity 

1.735ms  1 

Mass  fraction  of  H 2 

0.727 

Mass  fraction  of  H2O 

0.273 

Mixture  viscosity 

1.161  x  10~4kgm_1  s_1 

Re 

4 

Cathode 

Inlet  velocity 

7.33ms  1 

Mass  fraction  of  O2 

0.225 

Mass  fraction  of  N2 

0.751 

Mass  fraction  of  H2O 

0.024 

Mixture  viscosity 

2.87  x  10-5  kgm-1  s~l 

Re 

258 

Low  humidity 

Anode 

Inlet  velocity 

1.83  ms-1 

Mass  fraction  of  H2 

0.635 

Mass  fraction  of  H2O 

0.365 

Mixture  viscosity 

1.52  x  10~4kgm_1  s-1 

Re 

36 

Cathode 

Inlet  velocity 

7.91ms  1 

Mass  fraction  of  O2 

0.225 

Mass  fraction  of  N2 

0.734 

Mass  fraction  of  H2O 

0.046 

Mixture  viscosity 

3.71  x  10-5  kgm-1  s-1 

Re 

214 

High  humidity 

Anode 

Inlet  velocity 

2.21ms  1 

Mass  fraction  of  H2 

0.406 

Mass  fraction  of  H2O 

0.594 

Mixture  viscosity 

2.415  x  10~4 kgm-1  s-1 

Re 

3 

Cathode 

Inlet  velocity 

9.05  ms  1 

Mass  fraction  of  O2 

0.21 

Mass  fraction  of  N2 

0.705 

Mass  fraction  of  H2O 

0.085 

Mixture  viscosity 

5.194  x  10“5  kgm-1  s_1 

Re 

176 

Very  high  humidity 

Anode 

Inlet  velocity 

2.56  ms  1 

Mass  fraction  of  H2 

0.295 

Mass  fraction  of  H2O 

0.705 

Mixture  viscosity 

2.848  x  10~4kgm_1  s-1 

Re 

3 

Table  2  ( Continued ) 


Cathode 

Inlet  velocity 

12.9  ms-1 

Mass  fraction  of  O2 

0.187 

Mass  fraction  of  N2 

0.61 

Mass  fraction  of  H2O 

0.203 

Mixture  viscosity 

9.686  x  10-5  kgm-1  s-1 

Re 

134 

Temperature  of  cell 

343  K 

Pressure  of  cell 

1.013  x  105  Nm-2 

Critical  temperature  of  hydrogen  (7frH2 ) 

33.3  K 

Critical  temperature  of  air  (7"cr,air) 

132  K 

Critical  temperature  of  oxygen  (Tcr  o2) 

154.4  K 

Critical  temperature  of  water  (Tcr;Hoo) 

647.3  K 

Critical  pressure  of  hydrogen  (PCr,H2) 

12.8  atm 

Critical  pressure  of  air  (E’er, air) 

36.4  atm 

Critical  pressure  of  oxygen  (PCr,Oo) 

49.7  atm 

Critical  pressure  of  water  (PCr,Hoo) 

221.2  atm 

Io 

100.0  Am-2 

Dry  equivalent  weight  of  PEM  (Mm?dry) 

1.1kg  mol-1 

Open  circuit  potential 

1.1  V 

Cell  voltage 

0.53  V 

4.  Local  validation 

Comparison  of  the  predictions  for  the  present  model  with 
data  obtained  from  the  literature  is  only  a  partial  validation 
since  only  global  data  (e.g.  average  current  density)  could 
be  checked  against  experimental  measurements.  This  is  in¬ 
sufficient  as  it  does  not  give  enough  evidence  that  whatever 
predicted  in  the  modeled  fuel  cell  is  locally  correct.  Here,  val¬ 
idation  of  the  present  model  against  experiment  data  obtained 
using  a  segmented  cell  (to  allow  local  measurements  along 
the  channel)  will  be  carried  out.  The  segmented  cell  consists 
of  shoulder  widths  that  allow  data  such  as  the  current  density 
to  be  measured  locally  at  points  along  the  shoulders.  Details 
of  the  experimental  setup  can  be  found  in  the  work  of  Potter 
[36]. 

The  experiments  carried  out  used  hydrated  hydrogen  on 
the  anode  side  and  hydrated  air  on  the  cathode  side.  Experi¬ 
mental  data  were  obtained  for  different  cases  by  varying  the 
inlet  relative  humidity  and  stoichiometry  on  the  cathode  side 
(i.e.  relative  humidity  of  47%  at  stoichiometry  7,  relative  hu¬ 
midity  of  65%  at  stoichiometry  3,  relative  humidity  of  83% 
at  stoichiometry  2,  in  both  coflow  and  counterflow  configura¬ 
tions)  while  keeping  all  other  parameters  constant.  The  base 
case  conditions  used  for  the  experiments  are  summarised  in 
Table  3  [36]. 

A  cross-sectional  schematic  diagram  of  part  of  the  com¬ 
plete  fuel  cell  is  shown  in  Fig.  4.  The  repeating  unit  identi¬ 
fied  is  given  in  Fig.  5.  This  repeating  unit  will  be  the  solution 
domain  applied  for  the  present  model.  For  the  current  pre¬ 
diction,  the  following  grid  density  was  applied,  after  a  grid 
dependence  test  was  carried  out:  5  cells  in  the  /-direction  (x), 
10  cells  in  the j-direction  (y)  and  200  cells  in  the  k-direction  (z) 
for  the  half  size  gas  channels;  10  cells  in  the  /-direction  (x),  10 
cells  in  the  y-direction  (y)  and  200  cells  in  the  k-direction  (z) 
for  the  full  size  gas  channels;  and  49  cells  in  the  /-direction  (x), 
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Table  3 

Parameters  for  local  validation  of  model 


Channel  width 

1.0  x  10  3  m 

Channel  height 

1.0  x  10-3  m 

Channel  length 

0.137m 

Electrode  width 

6.0  x  10-3  m 

Electrode  height 

2.2  x  10-4  m 

Electrode  length 

0.137m 

Membrane  thickness  (tm) 

6.0  x  10-5  m 

Permeability  of  electrode  (at) 

1.0  X  10“ 10  m2 

Porosity  of  electrode  (e) 

0.4 

Anode  flowrate 

6.84  m/s 

Anode  mixture  density 

0.377  kg/m3 

Anode  mixture  viscosity 

2.106  x  10_4kgm_1  s-1 

Re  a 

12 

Case  1 

Cathode  stoichiometry  ( Rec  =  340) 

7 

%RHanode 

129 

%RHcalhodc 

47 

Cathode  mixture  density 

2.989  kgm-3 

Cathode  mixture  viscosity 

6.24  x  10-5  gm_1  s-1 

Cell  voltage 

0.5  V 

Total  current 

4.15  A 

Case  2 

Cathode  stoichiometry  ( Rec  =  253) 

3 

%RHanode 

129 

%RHcalhodc 

65 

Cathode  mixture  density 

3.115  kgm-3 

Cathode  mixture  viscosity 

3.16  x  10-5  kgm-1  s-1 

Cell  voltage 

0.5  V 

Total  current 

3.5  A 

Case  3 

Cathode  stoichiometry  ( Rec  =  132) 

2 

%RHanode 

129 

%RHcalh0dc 

83 

Cathode  mixture  density 

3.1  kg  m-3 

Cathode  mixture  viscosity 

3.548  x  10-5  kgm-1  s-1 

Cell  voltage 

0.5  V 

Total  current 

3.1  A 

Mole  fraction  ratio  of  O2/N2 

0.21/0.79 

Temperature 

333  K 

Pressure 

3.039  x  105 

P  sat 

20178. 96Nm-2 

Critical  temperature  of  hydrogen  (7fr 5h0  ) 

33.3  K 

Critical  temperature  of  air  (7"cr,air) 

132  K 

Critical  temperature  of  oxygen  (Tcr  o2) 

154.4  K 

Critical  temperature  of  water  (Tcr  Hio) 

647.3  K 

Critical  pressure  of  hydrogen  (Pcr.H? ) 

12.8  atm 

Critical  pressure  of  air  (PCr,air) 

36.4  atm 

Critical  pressure  of  oxygen  (PCr,Oo ) 

49.7  atm 

Critical  pressure  of  water  (Pcr.Hoo) 

221.2  atm 

to 

100.0  A  m“2 

Dry  density  of  PEM  (pm,dry) 

2000.0  kg  m-3 

Dry  equivalent  weight  of  PEM  (Mm^ry) 

1.1  kg  mol-1 

Open  circuit  potential 

1.1  V 

10  cells  in  the  j-direction  (y)  and  200  cells  in  the  ^-direction 
(z)  for  the  electrodes.  The  boundary  conditions  are  given  in 
Appendix  A. 

Experimental  data  has  been  taken  from  [36]  for  each  of 
the  cases  to  be  studied.  Since  the  value  of  porosity  is  not 
known  for  the  electrodes  used  in  the  experiments,  two  values 


of  porosity  are  presented  as  comparison.  The  current  den¬ 
sity,  species  mole  fractions  and  water  transport  coefficients 
presented  in  the  following  figures  are  obtained  at  the  centre 
of  the  cell  (i.e.  at  v  =  3.0  x  10-3  m,  y  =  1.22  x  10-3  m  at  the 
cathode). 

4.1.  Current  density  and  species  distribution 

Figs.  6  (a  and  b),  7(a  and  b)  and  8(a  and  b)  show  the  cur¬ 
rent  distribution  along  the  channel  for  a  cathode  inlet  relative 
humidity  of  47%,  65%  and  83%,  respectively,  for  both  the 
coflow  and  counterflow  configurations,  with  a  porosity  of  0.3 
and  0.4.  It  can  be  seen  that  a  porosity  of  0.3  gives  a  better  fit 
to  the  experiment  data  for  the  cases  of  65%  and  83%  relative 
humidity  at  the  cathode  while  it  under-estimates  the  current 
density  for  the  47%  relative  humidity  case.  This  is  true  for 
both  the  coflow  and  counterflow  configurations. 

With  a  porosity  of  0.4,  the  current  density  obtained  for 
the  47%  relative  humidity  case  shows  excellent  agreement  to 
that  obtained  from  the  experiment.  However,  they  are  higher 
for  the  cases  of  65%  and  83%  relative  humidity.  This  can  be 
explained  by  looking  at  the  mole  fraction  of  water  on  both 
the  anode  and  cathode  side.  From  Fig.  9(a-c),  it  can  be  seen 
that  water  concentration  at  the  anode  is  slightly  above  that 
for  the  liquid  water  saturation  value.  This  means  that  water  is 
present  in  the  liquid  phase.  On  the  cathode,  water  is  present 
in  the  liquid  phase  almost  along  the  whole  channel  length, 
except  near  the  inlet  where  water  is  present  in  the  vapour 
phase.  It  is  therefore  likely  that  liquid  water  is  reducing  the 
available  pore  volume  for  the  gas  flow,  especially  so  at  the 
cathode,  resulting  in  a  poorer  performance  for  all  three  cases 
considered. 

It  can  also  be  seen  that  at  a  higher  inlet  relative  humidity, 
a  higher  water  concentration  in  the  liquid  phase  at  the  cath¬ 
ode  is  present  along  the  whole  channel  length.  Since  in  the 
model,  water  is  only  considered  in  the  gas  phase,  phenom¬ 
ena  such  as  electrode  flooding  cannot  be  captured.  This  will 
result  in  a  higher  current  density  obtained  for  the  simulated 
results  compared  to  the  experimental  data  and  this  effect  will 
be  greater  as  the  relative  humidity  increases.  This  explains 
why  the  simulated  results  for  the  two  higher  humidity  cases 
produced  a  current  density  curve  which  is  higher  than  that 
obtained  from  the  experiment  while  the  simulated  results  for 
the  lowest  relative  humidity  case  shows  good  agreement  to 
the  experimental  data  when  a  porosity  of  0.4  is  used.  An¬ 
other  point  to  note  is  that  the  reduction  in  porosity  for  the 
two  higher  humidity  cases  can  effectively  represent  the  pore 
blockage  by  the  presence  of  liquid  water.  Hence,  a  porosity 
of  0.4  seems  to  be  a  better  fit  to  the  experimental  data  for  the 
47%  relative  humidity  case  while  a  porosity  of  0.3  is  a  better 
fit  for  the  two  higher  relative  humidity  cases. 

One  explanation  for  the  disagreement  near  the  outlet  may 
be  due  to  the  presence  of  liquid  water  in  the  electrodes.  Fiq- 
uid  water  tends  to  be  present  in  greater  amount  near  the  exit 
of  the  channel  since  water  is  accumulated  down  the  chan¬ 
nel  from  the  product  water  as  a  result  of  the  electrochemical 
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Fig.  4.  Cross-section  of  the  cell  geometry  used  in  the  experiments. 


reaction.  This  can  be  seen  in  [41-44]  which  show  that  the 
presence  of  liquid  water  does  lead  to  an  effective  pore  block¬ 
age.  The  lower  current  density  obtained  in  the  experiment 
may  be  explained  by  this  accumulation  of  water  that  caused 
mass  transport  limitation  of  the  reactant  gases  to  the  cata¬ 
lyst  sites  for  electrochemical  reaction  to  take  place.  Since  in 
the  model,  the  effect  of  two-phase  flows  has  not  been  im¬ 
plemented  and  all  the  water  produced  is  assumed  to  be  in 
the  vapour  state,  which  will  not  result  in  any  pore  blockage, 
such  two-phase  effects  will  not  be  visible.  Other  reasons  for 
the  disagreement  could  be  due  to  the  conditions  in  which  the 
experiments  have  been  conducted,  such  as  the  time  the  ex¬ 
periment  has  been  allowed  to  run  before  data  are  taken  and 
whether  the  MEA  has  been  used  a  couple  of  times  at  other  op¬ 
erating  conditions.  All  these  will  affect  the  experimental  data 
obtained. 

Referring  back  to  Fig.  9(a-c),  which  shows  the  distri¬ 
bution  of  species  along  the  channel  from  the  model,  it  can 
be  seen  that  due  to  the  high  stoichiometry,  the  mole  frac¬ 
tions  of  oxygen  and  hydrogen  remain  virtually  constant  in 
the  streamwise  direction  of  the  cell.  The  local  current  den¬ 


sity,  however,  varies  significantly  (e.g.  from  1.2  x  104Am  2 
to  8.8  x  103  Am  2).  Only  the  species  mole  fractions  in  the 
coflow  configuration  for  the  case  of  0.4  porosity  is  presented 
for  the  sake  of  brevity. 

4.2.  Cell  performance 

Returning  to  the  three  current  density  plots 
(Figs.  6  (a  and  b),  7(a  and  b)  and  8(a  and  b)),  the  case 
with  the  highest  stoichiometry  produced  a  higher  current 
density  along  the  whole  of  the  channel.  This  is  true  for  both 
the  coflow  and  counterflow  configurations.  The  average 
current  densities  obtained  along  the  channel  for  the  different 
cases  are  given  below  (at  a  cell  voltage  of  0.5  V): 

1.  At  a  relative  humidity  of  47%  (stoichiometry  =  7), 
the  average  current  density  is  1.0  x  104Am  2  and 
1.03  x  104  Am  2  for  the  coflow  and  counterflow  case, 
respectively. 

2.  At  a  relative  humidity  of  65%  (stoichiometry  =  3), 
the  average  current  density  is  8.0xl03Am~2  and 
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Fig.  5.  Repeating  unit  of  the  cell  used  in  the  experiment  and  model. 
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Fig.  8.  Local  current  density  down  a  cathode  side  gas  channel  at  a  cathode  relative  humidity  of  83%  in  (a)  coflow;  (b)  counterflow. 


9.0  x  103  Am-2  for  the  coflow  and  counterflow  case,  re¬ 
spectively. 

3.  At  a  relative  humidity  of  83%  (stoichiometry  =  2), 
the  average  current  density  is  7.0xl03Am~2  and 
8.0  x  103  Am-2  for  the  coflow  and  counterflow  case,  re¬ 
spectively. 

The  reason  for  this  is  because  a  higher  inlet  velocity  will 
allow  a  higher  concentration  of  reactant  gas  to  reach  the  cat¬ 
alyst  layer  for  electrochemical  reaction  to  take  place  by  im¬ 
proved  convection.  Furthermore,  with  a  lower  concentration 
of  water  fed  into  the  cell,  the  degree  of  mass  transport  limita¬ 
tion  occurring  in  the  electrodes  due  to  pore  volume  reduction 
by  liquid  water  will  be  reduced.  In  addition,  the  high  velocity 
helps  in  removing  any  excess  liquid  water  from  the  cell.  All 
these  factors  work  towards  a  better  performance  from  the  fuel 
cell,  provided  the  amount  of  water  introduced  into  the  cell  is 
still  sufficient  to  keep  the  membrane  well  hydrated.  Under  the 
present  operating  conditions,  the  counterflow  configuration 
tends  to  produce  a  marginally  higher  average  current  density 
compared  to  the  coflow  configuration. 

4.3.  Water  transport  in  the  membrane 

Fig.  10  shows  the  interaction  between  the  two  water  trans¬ 
port  mechanisms  in  the  membrane  of  a  fuel  cell:  electro- 
osmotic  drag  coefficient  and  diffusion  coefficient  at  a  cath¬ 
ode  relative  humidity  of  47%,  65%  and  83%.  The  net  trans¬ 
port  of  water  per  molecule  of  proton  across  the  membrane 


is  also  shown.  The  net  transport  of  water  per  proton  is  al¬ 
ways  a  positive  number.  This  means  that  there  is  a  positive 
flow  of  water  from  the  anode  to  the  cathode  along  the  whole 
channel.  In  addition,  overall  diffusion  takes  place  from  the 
cathode  to  the  anode  along  the  whole  channel.  At  the  inlet, 
the  amount  of  reactants  and  the  partial  pressure  of  water  are 
high,  leading  to  high  membrane  conductivity  and  hence  the 
local  current  density  is  also  high.  Further  down  the  channel, 
the  membrane  becomes  dryer  due  to  lower  water  concen¬ 
tration  on  the  anode  side.  This  causes  the  membrane  con¬ 
ductivity  and  local  current  density  to  decrease.  Decreases  in 
local  current  density  then  lead  to  lower  water  concentration 
on  the  cathode  side  and  hence  the  diffusion  coefficient  of 
water  from  the  cathode  to  anode  decreases.  The  net  water 
transport  per  proton  decreases  down  the  channel  since  there 
is  less  water  present  down  the  channel  and  hence  the  amount 
that  each  proton  can  drag  with  it  across  the  membrane  de¬ 
creases. 


5.  Parametric  studies 

The  following  sections  will  discuss  the  performance  of  a 
PEMFC  upon  variation  in  geometrical  and  operating  condi¬ 
tions,  based  on  Fig.  2.  To  elucidate  some  of  the  major  con¬ 
tributions/limitations  (i.e.  cell  design  and  material  proper¬ 
ties),  the  following  cases  were  studied:  the  effect  of  elec¬ 
trode  thickness;  the  effect  of  shoulder  width  with  respect  to 
the  gas  channel  width;  the  permeability  effect  of  the  elec- 
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trodes;  the  effect  of  using  oxygen  instead  of  air.  The  base 
case  conditions  to  be  used,  which  are  typical  for  fuel  cell  op¬ 
erations,  are  given  in  Table  4.  Otherwise,  those  from  Table  3 
are  used.  Only  one  parameter  was  changed  from  the  base  case 
conditions  at  a  time.  The  boundary  conditions  are  given  in 
Appendix  A. 

5.7.  Parametric  variations  relevant  to  cell  geometry 

In  many  fuel  cell  applications,  it  is  important  to  consider 
the  geometry  of  the  cell  so  that  maximum  current  can  be 
obtained  from  as  small  a  cell  as  possible.  This  will  reduce 
the  cell  size,  especially  in  fuel  cell  stacks,  which  will  make 
the  design  more  attractive.  Hence,  considerations  have  to  be 
made  in  terms  of  how  many  gas  channels  there  should  be  with 
respect  to  the  shoulder  width  in  a  given  area  of  cell  space, 
and  the  thickness  of  the  electrodes  to  be  used.  These  factors 
will  be  considered  below. 


5.7.7.  Effect  of  shoulder  width 

To  study  this  effect,  the  gas  channel  width  is  kept  at 
1.0  x  10-3  m  while  the  width  of  the  shoulder  is  varied  be¬ 
tween  7.5  x  10  4  m  to  2.5  x  10  4  m.  The  size  of  the  shoul¬ 
der  width  affects  the  performance  of  the  fuel  cell  in  three 
ways.  Firstly,  having  a  smaller  shoulder  width  will  enhance 
the  transport  of  reactant  species  to  the  catalyst  layer  directly 
above/below  the  shoulder  area.  This  can  be  clearly  seen  in 
Fig.  11.  Fig.  1 1  shows  the  distribution  of  oxygen  at  the  cath¬ 
ode  electrode/membrane  interface  when  the  shoulder  width 
is  varied  while  keeping  the  channel  width  constant.  As  the 
shoulder  width  increases,  the  oxygen  concentration  over  the 
shoulder  area  decreases  significantly.  This  tends  to  decrease 
the  limiting  current  density.  Furthermore,  the  oxygen  con¬ 
centration  over  the  channel  area  is  lower  since  some  of  the 
oxygen  tends  to  diffuse  to  the  wider  shoulder  area. 

Secondly,  having  a  smaller  shoulder  width  will  increase 
the  contact  resistance  between  the  current  collector  and  the 


K.W.  Lum,  J.J.  McGuirk  /  Journal  of  Power  Sources  143  (2005)  103-124 


115 


c 

o 


(a)  o  Channel  length,  m 


c 

CD 

O 

a> 

o 

o 

c 
o 
c n 

3 


Q 


c 

o 


c 


Channel  length,  m 


Fig.  10.  Electro-osmotic  drag  coefficient,  diffusion  coefficient  and  net  water  transport  across  the  membrane  for  a  cathode  relative  humidity  of  (a)  47%;  (b) 
65%;  (c)  83%  in  coflow. 


ME  A.  This  will  lead  to  greater  ohmic  loss  since  there  is 
less  contact  area  to  collect  the  electrons.  However,  the  cur¬ 
rent  density  in  Fig.  12  shows  that  as  the  shoulder  width 
increases,  the  average  current  density  obtained  is  lowered 
(from  1.27  x  104  Am-2  to  1.1  x  104  Am-2)  as  a  result  of 
the  highly  reduced  oxygen  concentration  above  the  shoulder 
area.  Hence,  for  the  present  operating  conditions,  it  is  advan¬ 
tageous  to  have  a  higher  oxygen  concentration  at  the  cathode 
electrode/membrane  interface  by  having  a  smaller  shoulder 
width,  provided  the  properties  of  the  porous  backing  outlined 
below  are  not  adversely  affected.  Thirdly,  the  ratio  between 
shoulder  and  channel  widths  will  affect  the  clamping  pres¬ 
sure  distribution  that  the  porous  backing  experiences,  which 


in  turn  affects  properties  such  as  porosity  distribution,  ther¬ 
mal  impedance  and  contact  resistances  [45]. 

The  polarization  curves  shown  in  Fig.  1 3(a)  show  that  with 
a  wider  shoulder  width,  mass  transport  limitation  takes  place 

Table  4 


Parameters  for  parametric  studies 


Channel  length 

0.1  m 

Electrode  width 

1.6  x  10~3  m 

Electrode  height 

2.0  x  10~4  m 

Porosity  of  electrode  (e) 

0.3 

%RHcalh0dc 

83 

Cathode  mixture  viscosity 

3.36  x  10~5  kgm-1  s_1 

Total  current 

3.1  A 
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at  a  lower  current  density,  resulting  in  a  lower  maximum 
peak  for  the  power  density  curves,  Fig.  13(b).  Having  a  ratio 
of  gas  channel  width  to  shoulder  width  of  1.3  (for  the  case  of 
shoulder  width  equivalent  to  7.5  x  10-4  m)  is  seen  to  produce 
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Fig.  12.  Average  current  density  at  the  electrode/membrane  interface  for  the 
different  shoulder  widths. 


the  lowest  peak  performance.  Increasing  the  ratio  to  2  shows 
an  increase  in  the  peak  performance.  Further  increment  of  the 
ratio  to  4  results  in  a  performance  only  slightly  better  than 
that  using  a  ratio  of  2.  Hence,  optimum  gas  channel  width 
to  shoulder  width  is  seen  to  occur  at  a  ratio  of  approxima¬ 
tely  2. 


5.7.2.  Effect  of  electrode  thickness 

Four  cases  of  electrode  thicknesses  are  considered, 
namely:  1.0xl0_4m,  2.0xl0_4m,  4.0xl0_4m  and 
6.0  x  10~4  m.  As  the  electrode  thickness  increases,  a  lower 
concentration  of  oxygen  is  present  at  the  reactive  area,  as  can 
be  seen  in  Fig.  14.  This  might  easily  lead  to  oxygen  depletion 
at  high  current  density,  especially  near  the  end  of  the  channel. 
However,  for  a  thicker  electrode,  there  is  a  more  even  distri¬ 
bution  of  oxygen  concentration  between  the  channel  area  and 
the  shoulder  area  since  there  is  more  room  for  diffusion  in 
the  span  wise  direction. 

Fig.  15  shows  the  average  current  density  obtained  down 
the  channel  for  the  different  electrode  thicknesses.  Near  the 
inlet,  the  thinner  electrode  produces  higher  current  density. 
This  is  because  a  higher  concentration  of  oxygen  reaches  the 
reaction  area.  However,  the  thicker  electrode  outperforms  the 
thinner  one  10%  down  the  channel.  This  is  clearly  not  due  to 
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Fig.  13.  (a)  Polarization  curve;  (b)  power  density  curve  for  the  various  shoulder  widths. 
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a  reduction  in  the  amount  of  oxygen  in  the  thinner  electrode. 
The  reason  may  be  seen  in  Fig.  16. 

Fig.  16  shows  the  amount  of  water  present  at  the  anode 
electrode/membrane  interface.  As  can  be  seen,  about  10% 
down  the  channel  from  the  inlet,  there  is  less  water  avail¬ 
able  on  the  thinner  electrodes.  This  results  in  membrane  de¬ 
hydration  and  hence  reducing  the  membrane  conductivity. 
As  a  result,  the  electrokinetics  is  limited  by  the  transport 
of  protons  across  the  membrane  to  the  cathode,  leading  to 
lower  current  density.  This  has  happened  despite  the  fact  that 
there  is  sufficient  reactant  gas  at  both  the  anode  and  cathode 
for  electrochemical  reaction  to  take  place  (Figs.  14  and  17). 
From  these  results,  the  importance  of  keeping  the  membrane 
well  hydrated  to  sustain  high  performance  is  clearly  shown. 
Increasing  the  electrode  thickness  is  seen  to  improve  the 
cell  performance  to  the  point  where  mass  transport  limita¬ 
tion  begins  to  take  place.  This  occurs  for  the  case  where  the 
electrode  thickness  reaches  6.0  x  10-4  m.  There,  despite  the 
fact  that  water  is  present  in  higher  amounts,  hydrogen  does 
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Fig.  16.  Water  distribution  at  the  anode  electrode/membrane  interface  at  different  electrode  thicknesses  of:  (a)  1 .0  x  10  4  m;(b)  2.0x10  4  m;(c)  4.0x10  4  m; 
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not  reach  the  catalyst  layer  in  an  amount  required  for  effi¬ 
cient  electrochemical  reaction  to  take  place,  as  can  be  seen  in 
Fig.  17. 

Fig.  18(a  and  b)  show  the  polarization  and  power  den¬ 
sity  curves  for  the  cell  with  different  electrode  thick¬ 
nesses  respectively.  Increasing  the  electrode  thickness  from 
1.0  x  10  4  m  to  2.0  x  10  4  m  results  in  better  performance 
(average  current  density  increases  from  1.17  x  104Am  2 
to  1.2  x  104Am-2).  Improvement  is  seen  on  further  in¬ 
creasing  the  thickness  to  4.0  x  10-4  m  (average  current  den¬ 
sity  =  1 .22  x  104  A  m~2).  Mass  transport  limitation  is  seen  to 
take  place  when  a  thickness  of  6.0  x  10-4  m  is  used,  where 
the  polarization  curve  drops  below  the  level  in  which  a  thinner 
electrode  thickness  was  used.  In  general,  a  thinner  electrode 
thickness  is  preferable  as  it  improves  the  fluxes  of  reactant 
gases  to  the  catalyst  layer.  The  problem  of  dehydration  at  the 
anode  can  be  easily  reduced  by  sufficiently  humidifying  the 
inlet  fuel  so  that  the  membrane  can  be  kept  hydrated  to  pre¬ 
vent  drying.  In  this  case,  an  electrode  thickness  of  between 
2.0  x  10-4  m  to  4.0  x  10-4  m  is  preferable. 


5.2.  Parametric  variations  relevant  to  operating 
conditions 

Apart  from  geometrical  considerations,  different  operat¬ 
ing  conditions  will  also  affect  the  performance  of  the  fuel 
cell.  Factors  such  as  electrode  permeability,  operation  using 
air  or  oxygen,  among  others,  are  believed  to  play  a  part.  The 
role  of  all  these  factors  can  easily  be  studied  using  the  cur¬ 
rent  model.  In  this  section,  an  analysis  of  such  parameters 
has  been  carried  out. 

5.2.7.  Effect  of  electrode  permeability 

Three  different  values  of  permeability  are  considered  here: 
1.0  x  10~8  m2,  1.0  x  10~9  m2  and  1.0  x  10~10  m2.  The  role 
of  the  degree  of  permeability  only  affects  the  flux  of  species 
towards  the  catalyst  layer  by  convection,  due  to  the  Darcy 
term.  As  the  permeability  value  decreases,  there  is  more  re¬ 
sistance  to  the  flow  of  species  to  the  catalyst  layer.  This  can 
be  seen  in  Fig.  19.  Decreasing  the  electrode  permeability  al¬ 
lows  the  reactant  species  to  spread  over  to  the  shoulder  area. 
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(b)  Current  density,  A/m2 

Fig.  18.  (a)  Polarization  curve;  (b)  power  density  curve  for  the  various  electrode  thicknesses. 
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Fig.  20.  Average  current  density  at  the  electrode/membrane  interface  for  the 
different  permeability  values. 


This  can  enhance  the  gases  to  flow  in  the  spanwise  direction 
to  better  utilise  the  catalyst  layer  above  the  shoulder  areas. 
However,  this  will  only  improve  the  cell  performance  if  the 
shoulder  area  is  greater  than  the  channel  area.  As  mentioned 
above,  the  use  of  a  gas  channel  having  a  smaller  width  com¬ 
pare  to  the  shoulder  is  not  preferable.  The  abrupt  changes  that 
can  be  discerned  in  Fig.  19(a)  close  to  the  outlet  originate 
from  the  very  high  absolute  permeability  (1.0  x  10  8m2) 
in  the  porous  backing,  leading  to  a  higher  convective  flow 
in  the  streamwise  and  spanwise  directions  than  for  the  two 
cases  with  lower  permeability.  Furthermore,  the  use  of  the 
outlet  boundary  condition  here  requires  the  flow  to  be  fully 
developed,  which  is  not  the  case  here  for  such  a  high  value  of 
permeability.  However,  since  the  error  is  small,  it  is  admis¬ 
sible  to  retain  this  boundary  condition  as  the  abrupt  changes 
are  only  visible  for  the  high  permeability  case. 


The  corresponding  average  current  density  for  the  three 
cases  is  given  in  Fig.  20.  It  can  be  seen  that  between  the  per¬ 
meability  value  of  1.0  x  10-8  Am-2  to  1.0  x  10  9  Am  2,  a 
drop  in  average  current  density  is  obtained  (1.28  x  104Am  2 
to  1.205  x  104  Am  2).  As  the  permeability  value  decreases 
further,  no  significant  decrease  in  average  current  density 
is  observed  (1.202  x  104  Am-2).  Hence,  there  is  a  “limit¬ 
ing”  permeability  value  whereby  further  reduction  will  not 
affect  the  overall  performance  of  the  fuel  cell.  This  is  fur¬ 
ther  emphasised  in  the  polarization  and  power  density  curves 
in  Fig.  21  (a  and  b)  respectively.  In  general,  it  may  be  sug¬ 
gested  that  in  terms  of  electrode  properties,  the  porosity 
term  may  be  more  important  in  affecting  the  performance 
of  a  fuel  cell.  This  means  that  diffusion  plays  a  more  im¬ 
portant  role  than  convection  to  transport  the  specie  fluxes 
to  the  catalyst  layer.  From  this  study,  it  can  be  seen  that 
lowering  the  permeability  of  the  electrode  beyond  a  value 
of  1.0  x  10  9m2  will  not  affect  the  performance  of  the 
cell. 

5.2.2.  Effect  of  oxidant  concentration 

In  fuel  cell  operations,  it  is  common  to  use  air  as  the  ox¬ 
idant  as  air  is  readily  available  and  does  not  need  to  be  pro¬ 
cessed.  However,  the  concentration  of  oxygen  in  air,  which 
is  the  main  reactant  gas,  is  low.  Performance  will  be  affected 
unless  a  sufficiently  high  pressure  is  applied  at  the  cathode 
inlet  to  increase  the  concentration  of  oxygen  reaching  the 
catalyst  layer.  This  is  because  the  current  density  is  depen¬ 
dent  on  the  partial  pressure  of  oxygen.  Fig.  22  compares  the 
performance  of  the  cell  by  using  different  concentrations  of 
oxygen  as  the  oxidant.  The  three  cases  considered  are:  air; 
air  enriched  with  30%  of  oxygen;  pure  oxygen. 

It  can  be  seen  that  not  much  difference  is  observed  be¬ 
tween  all  the  three  cases  at  high  potential  since  low  current 
density  is  produced  and  the  need  for  efficient  distribution  of 
oxygen  to  the  catalyst  layer  has  not  set  in.  At  low  potential, 
mass  transport  limitation  begins  to  set  in  for  the  case  where 
pure  air  is  used  (^0.7  V).  This  is  due  to  insufficient  oxygen 


(a)  Current  density,  A/m2 


Fig.  21.  (a)  Polarization  curve  (b)  power  density  curve  for  the  various  permeability  cases. 
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Fig.  22.  Polarization  curves  for  the  different  concentration  of  oxygen  ap¬ 
plied. 

reaching  the  catalyst  layer  where  electrochemical  reaction 
takes  place.  With  oxygen-enriched  air,  mass  transport  limi¬ 
tation  takes  place  at  a  lower  potential  0.6  V)  and  hence 
higher  current  density  is  possible.  When  pure  oxygen  is  used 
as  the  oxidant,  the  polarization  curve  obtained  is  even  higher 
since  at  the  same  applied  pressure,  a  higher  concentration  of 
oxygen  reaches  the  catalyst  layer.  No  mass  transport  limita¬ 
tion  is  observed  for  the  case  with  pure  oxygen  as  the  oxidant. 
The  factor  that  prevents  the  cell  from  producing  higher  cur¬ 
rent  when  pure  oxygen  is  used  is  probably  due  to  the  ohmic 
losses  that  occur  at  the  membrane.  Hence,  the  use  of  pure 
oxygen  for  fuel  cell  application  is  preferable  although  it  has 
to  be  weighed  against  the  added  cost  in  doing  so. 

6.  Conclusions 

A  three-dimensional  model  of  a  complete  polymer  elec¬ 
trolyte  membrane  fuel  cell  was  implemented  in  an  in-house 
code,  which  has  been  validated  on  a  global  basis  with  the 
work  of  Shimpalee  et  al.  [22].  From  the  results  obtained,  the 
presence  of  liquid  water  on  the  cathode  side  is  seen  to  lead 
to  blockage  of  the  electrode,  causing  mass  transport  limita¬ 
tion  to  the  catalyst  layer  for  electrochemical  reaction  to  take 
place.  This  will  result  in  a  lower  current  density  obtained. 

The  model  was  then  validated  on  a  local  basis  with  ex¬ 
perimental  data  of  current  density,  which  were  obtained  lo¬ 
cally  along  the  channel  length.  This  is  a  stronger  validation 
as  it  verifies  current  density  obtained  point  by  point  along 
the  channel  length.  Excellent  agreement  has  been  obtained 
between  experimental  data  and  the  simulated  results  locally. 

The  net  transport  of  water  across  the  membrane  from  the 
anode  to  the  cathode  side  was  seen  to  decrease  down  the 


channel  length.  This  is  because  the  reactant  concentration  de¬ 
creases  down  the  channel,  hence  the  current  density  obtained 
will  be  lower,  leading  to  lower  water  produced  on  the  cath¬ 
ode.  This  will  reduce  the  effect  of  back  diffusion,  leading  to  a 
lower  concentration  of  water  on  the  anode  side  to  be  dragged 
by  the  protons  across  the  membrane.  Overall  transport  of  wa¬ 
ter  is  seen  to  take  place  from  the  anode  to  the  cathode.  As 
a  result,  it  is  crucial  to  humidify  the  anode  reactants  to  keep 
the  membrane  well  hydrated  to  prevent  dehydration. 

Parameter  variations  relevant  to  the  design  of  fuel  cells 
have  been  carried  out  to  study  their  effect  on  the  cell  oper¬ 
ation.  It  has  been  found  that  a  thinner  electrode  layer,  and  a 
smaller  shoulder  to  channel  width,  is  preferable  since  reactant 
gases  are  able  to  reach  the  catalyst  layer  with  less  resistance. 
Lowering  the  permeability  of  the  electrode  enhances  species 
transport  in  the  spanwise  direction  but  this  is  only  useful  if 
the  shoulder  area  is  large,  which  is  not  preferred  as  in  doing 
so  prevents  the  reactant  gases  from  being  efficiently  trans¬ 
ported  to  the  catalyst  layer.  Finally,  introducing  pure  oxygen 
as  the  oxidant  produces  much  higher  current  density  (at  the 
onset  of  mass  transport  limitation)  compared  to  when  pure 
air  or  air  enriched  with  oxygen  is  used  as  the  oxidant. 

In  general,  there  are  many  factors  affecting  the  perfor¬ 
mance  of  the  fuel  cell  simultaneously  and  all  these  factors 
have  to  be  weighed  comparatively  together,  depending  on 
the  kind  of  output  desired  for  a  particular  application.  It  has 
been  demonstrated  in  this  work  that  it  is  relatively  easy  and 
straightforward  to  use  the  current  CFD  model  when  design¬ 
ing  fuel  cells  to  test  various  parameters  to  obtain  the  optimum 
operating  conditions  and  cell  geometry.  A  plot  of  the  power 
density  curve  will  show  the  condition  under  which  maximum 
performance  occurs.  This  will  aid  in  designing  and  optimising 
fuel  cell  performance  without  the  need  to  carry  out  expensive 
and  time-consuming  experiments. 
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Appendix  A 

Boundary  conditions 

The  boundary  conditions  are  denoted  with  roman  numer¬ 
als  for  easy  identification  in  the  two  geometries  considered 
here,  as  shown  in  Figs.  2  and  5. 

A.l.  Inlet  (I) 

At  the  channel  inlet,  the  known  inlet  velocity  vectors  and 
species  mass  fractions  are  specified.  An  “in”  superscript  rep¬ 
resents  condition  at  the  inlet  of  the  channel  and  the  subscript 
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“i”  refers  to  the  individual  species, 
u  =  um,  mi  =  m)n 

A.2.  Outlet  (II) 

At  the  channel  outlet,  zero  gradient  is  applied  to  the  ve¬ 
locity  vectors  and  species  mass  fractions, 

du  dmi 

—  =  0,  — -  =  0 

dz  dz 


and  fluxes  of  oxygen  and  water  are  continuous  across  the 
interface, 


(m)p, 


dmf 
3 y 


ji  3(w)p 

€  dy 


du 


where  the  superscript  “np”  refers  to  properties  in  the  non- 
porous  region  while  p  refers  to  properties  in  the  porous  region. 


A3.  Wall  (III) 


For  a  porous  media  with  an  impermeable  surface,  the  ve¬ 
locity  normal  to  the  surface  is  zero  while  the  velocities  in  the 
other  directions  have  a  “slip”  condition  whereby  their  gradi¬ 
ent  is  zero.  Considering  such  a  circumstance  for  the  porous 
electrodes:  In  the  X-Z  plane, 


dmi 

dy 


In  the  X-Y  plane, 


For  the  non-porous  impermeable  surfaces,  considering 
boundaries  in  the  X-Z  plane, 


dmi 

dy 


In  the  X-Y  plane, 


dmi 

dx 


A.4.  Symmetry  (IV) 


Assuming  that  end  effects  are  negligible  and  that  only  a 
repeating  unit  of  the  complete  cell  needs  to  be  modeled,  each 
repeating  unit  will  be  bounded  by  two  symmetry  planes.  For 
the  porous  electrodes,  symmetries  occur  in  the  Y-Z  plane, 

d(u)  dmi 

—  =0,  — -  =  0 
dx  dx 

For  the  non-porous  channels, 


du 

dx 


dmi 

dx 


A.  5.  Porous/non-porous  interface  (V) 

For  complex  geometries,  different  blocks  have  to  be  placed 
together.  At  the  interface  between  the  porous  backing  and  the 
flow  channels,  the  point  wise  velocities  and  pressure  in  the 
plain  fluid  (flow  channels)  must  be  coupled  with  their  super¬ 
ficial  counterparts  in  the  porous  medium.  The  mass  fractions 


A.  6.  Membrane/catalyst  interface  (VI) 

The  electrochemistry  at  the  membrane/catalyst  interfaces 
on  the  cathode  and  anode  side  is  implemented  via  source 
terms  for  the  computational  nodes  adjacent  to  the  membrane, 
as  given  in  Table  1 . 
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